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ABSTRACT 

This paper presents applications of weighted meshless scheme for conservation laws 
to the Euler equations and the equations of ideal magnetohydrodynamics. The diver- 
gence constraint of the latter is maintained to the truncation error by a new meshless 
divergence cleaning procedure. The physics of the interaction between the particles is 
described by an one-dimensional Riemann problem in a moving frame. As a result, 
necessary diffusion which is required to treat dissipative processes is added automat- 
ically. As a result, our scheme has no free parameters that controls the physics of 
inter-particle interaction, with the exception of the number of the interacting neigh- 
bours which control the resolution and accuracy. The resulting equations have the 
form similar to SPH equations, and therefore existing SPH codes can be used to im- 
plement the weighed particle scheme. The scheme is validated in several hydrodynamic 
and MHD test cases. In particular, we demonstrate for the first time the ability of a 
meshless MHD scheme to model magneto-rotational instability in accretion disks. 
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1 INTRODUCTION 

Computational magnetohydrodynamics (MHD) remains an important tool to understand complex behaviour of astrophys- 
ical plasmas. While many methods have been developed to solve equations of ideal MHD of Eulerian (cartesian) meshes, 
few successful Lagrangian meshless formulations exists. The latter, however, are desired for problems which lack particular 
symmetries, cover many length-scales or require adaptivity, for example stellar collisions star or star cluster formations. 

Smoothed particle hydrodynamics (SPH) proved to be a successful Lagrangian meshless scheme to solve equations of fluid 
dynamics in variety research flelds ( Monaghan|2005" I . Despite its limitations, it has been also successfully used in wide range 
of astrophysical problems, including evolution of gaseous disks around black holes or stars, star formation, stellar collisions, 
and cosmology. Because of its simplicity and versatility, several attempts, albeit with limited success, have been made to 



include magnetic flelds into SPH, thereby formulating smoothed particle MHD, or SPMHD for short (Price & Monaghan 
|2OO5[|B0rve et al.|2006[ [Rosswog fc Price|2007| . It soon became clear that SPMHD equations where plagued with two main 
problems: tensile instability and maintenance of the divergence constraint, V ■ B = 0. The former is a general problem of SPH 



namely the equations are unstable to tensile stresses ( Swegle et al.|1995 Monaghan|2000 1 . In case of SPMHD, this instability 
manifest itself as particle clumping in the regions where magnetic pressure dominates gas pressure, and therefore rendering the 
simulation of strongly magnetised plasma unfeasible. It has been shown that the equations can be stabilised by either adding 
short-range repulsive forces between particles (e.g. Price Sz Monaghan|2004 1 , or sacriflcing momentum conservation (e.g. Price 
& Monaghan 2005). Alternatively, B0rve et al. (2001 1 showed that the stability of SPMHD equations can also be achieved by 



adding a source term proportional to the divergence of magnetic fleld to the momentum equation. While these approaches 
appear to remove the tensile instability, the divergence constraint still remains an issue in these SPMHD formulations. [Rosswog| 
fc Price ( 2007 1 were able to formulate manifestly divergence-free SPMHD equations which appears to work for a wide range 
of problems (e.g. [Price fc Rosswog|[2^006[ [Price fc Bate||2008[ ). However, they solve a limited form of the induction equation 
which only permits topologically trivial fleld conflgurations and is unable to model more complicated MHD phenomena. 



such as magneto-rotational instability. Dolag fc Stasyszyn ( 2009 1 were also able to formulate stable SPMHD equations by 



using a combination of several techniques, such as the addition of the source term proportional to the magnetic divergence 
to the momentum equation, artiflcial dissipation and smoothing of the magnetic field, and modification of the induction 
equation. This approach introduces several free parameters, which the authors were able to constrain by fitting their results 
to the solutions of several shock tube problems computed with conservative Eulerian MHD schemes. Despite the progress 
in this field, the difficulties associated with formulating consistent SPMHD equations advocate the need of the alternative 



approaches to formulate meshless Lagrangian MHD schemes. In a gradient particle magnetohydrodynamics (GPMHD, Maron 
[fc Howes[[2003[ [Maron[[2005[ ), the equations of ideal MHD were discretised on a set of particles by fltting a second or forth 
order polynomial into the data in order to obtain first and second order derivatives of the desired quantities. However, such 
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approach still requires addition of artificial diffusion to the resulting equations to model dissipative and resistive processes 
across discontinuities, and this introduces additional free parameters which control dissipative processes. In addition, the local 
conservation, for example third Newton law, is satisfied to the truncation error only. This is a potential source of error for 
the problems with interacting strong shock waves, since the truncation error across a shock wave is always of order of unity. 
Nevertheless, GPMIfD appears to be a viable, albeit noticeably more complex, alternative for SPMHD, at least for subsonic 
and weakly supersonic flows. 

Alternatively, one may attempt to utilise a Godunov approach. It has been shown that most of, if not all, SPH limitation 
can be removed by solving a hydrodynamic Riemann problem between each pair of interacting particles to obtain pressure 
forces, instead of computing separately pressure forces and artificial viscosity terms ( ^Inutsuka 2002j Cha fc Whitworth,20d3l 
|Chaet al.|20T0l ). In such Godunov SPH (GSPH) formulations, the necessary mixing and dissipation is included into underlying 
Riemann problem which described the physics of the interaction. Borrowing these ideas, it is therefore conceivable that 
Godunov-like meshless MHD formulation will eliminate the problems that plague SPMHD formulation, thereby permitting 
formulation of consistent Lagrangian meshless MHD scheme. 

In this paper we formulate a weighted particle MHD scheme. Our scheme is based on a meshless discretisation of the 
conservation law equations, which was pioneered by Vila (19991, and in Section [2] we present a heuristic derivation of the 
meshless conservative equations. We give the implementation details of our scheme in Section |3] In Section [4.1| we present 
applications to the equation of ideal hydrodynamics, and in Section |4.2| we show how our scheme can be applied to the 
equations of ideal MHD. In Section |5] we validate our meshless MHD scheme on several test problems, and finally, we present 
our conclusions in Section IH 



2 METHODS 

2.1 Meshless equations for conservation laws 

In what follows, we present a heuristic derivation of meshless discretisation of a scalar conservation law. The readers interested 
in a rigorous mathematical formulation supplemented with convergence theorems are referred to the original papers by | Vila| 
(19991 and |Lanson fc VHa ( 2008a|b' I. Following these works, a weak solution to a scalar conservation law 

|^+V-(F + an) = S, (1) 

is defined by 

J {u{x,t)ip + ¥{u,:>c,t) ■\/(p + S{:>i,t)ip) d:>idt = 0. (2) 

Here, m(x, t) is a scalar field, ^(x, t) is its source, F(ii, x, t) is its fiux in a frame moving with velocity a(x, t), and the integral is 
carried out over all space-time domain of a problem at hand. The function ip = f{x., t) is an arbitrary differentiable function in 
space and time, ip = dip/dt + a(a;, t) ■ Vtp is an advective derivative, and a(x, t) is an arbitrary smooth velocity field which will 
describe motion of particles. This integral is discretised on a set of particles with coordinates Xi with the help of a partition 
of unity 

V'.(x) = u-(x)H/(x-x,,/i(x)), (3) 

where, u'(x)^^ — IV(x — Xj , /i(x)) is an estimate of the particle number density and the sum is carried out over all particles, 
/i(x) is a smoothing length, and W['x., h) is a smoothing kerneQwith a compact support of size h; in what follows we assume 
that the kernel is normalised to unity. Inserting 1 = J]]; ''Pii'^) into an integral of an arbitrary function, we obtain 



(4) 



where Vi = J i/'i(x) dx is the effective volume of a particle i, and in the third term we use first-order Taylor expansion of /(x). 
In principle, a higher order discretisation is also possible, but for the purpose of this work such an one-point quadrature is 



sufficient; in fact, on a regular distribution of particles this discretisation is second order accurate (c.f. [3.21. Application of 
this discretisation to Eq. [2] gives 

f (V^u^if, + V,Fr{D°'p)^ + V,S^^,) = 0. (5) 

Here, the Einstein summation is assumed over Greek indexes, which refer to the components of a vector. Also, the gradient of 
a function Lp at i-particle location, {Vip)f is replaced by its discrete version {D°'ip)i, and this, for example, can be computed 
with SPH estimates of a gradient. However, a much better approach is to employ a more accurate meshless gradient estimate 



^ In this paper we use cubic-spline smoothing kernel, which is commonly used in SPH. 
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suggested by Lanson & Vila (2008a I. They showed that a second order accurate meshless partial derivative is given by the 
following expression 

{D-f)^ = ^(/, - /,)i3r''A4v.(xO = ^(/, - /OV'^xO, (6) 

J 3 

where V'f (xi) = -Bf Ax^-t/jj (xi), Aa;^ = (xj — Xi)'* and B"^ — {E'^^)~^ is a renormalisation matrix defined by its inverse 

i5r'' = ^A:r5A4v.(x.)- (7) 

j 

Finally, integrating the first term by parts with assumption that vanishes at boundaries, and applying the following rear- 
rangement to the second term 

permits separation of <^ from the rest 

j dtJ2^^ (^|(^""^) - E [V^Ft^n^^) - V,J;"C(X,)] + V.S.^ = 0. (9) 

The above is true for an arbitrary function ip if the expression in brackets vanishes, namely 

f^iv.u,) + J2 [y-^F"^^7i^^) - v,F;'ijn^3)] = (10) 

The extension of this equation to a general vector field u is straightforward: this equation is applied to each component of the 
field. These equations are similar to SPH equations, with the difference that the physics of the particle interaction is hidden in 
the fluxes and source terms. Indeed, if one uses the fluxes of ideal Lagrangian hydrodynamics, the equation of motions similar 
to SPH can be derived. As with SPH, such equations do not include dissipative processes, and therefore must be augmented 
with explicit diffusive terms, namely in the form of artificial viscosity, conductivity and resistivity. 

However, the power of the new scheme becomes apparent with the realisation that one can utilise the fluxes produced by 
the solution of an appropriate Riemann problem between particles i and j, which automatically include necessary dissipation. 
Defining such a flux as F^, and setting F"" — F" — F^ gives 

+ [K^^x.) - V,-C(x,)] = S,K. (11) 

3 

Finally, deflning a vector nfj = Vitpjlxi) — Vj-^f (x_,) and n = n/|n|, the equations take the following form 

^(Vi^O + ^(F.j ■ n.,)|n,,l = S^V^. (12) 

3 

It becomes clear, only the projection of the flux on the direction of the vector riij is required, and therefore for a wide range 
of problems the flux can be obtained by solving an appropriate ID Riemann problem in a frame moving with mean velocity 
of the two particles, a^j . 



2.2 Linear monotonic reconstruction 



The one dimensional flux in Eq. [T2]is naturally computed at the midpoint between particles i and j, i.e. at Xij — (x^ -|-Xj)/2. 
To achieve higher than the first order accuracy, it is necessary to linearly reconstruct left and right states of the Riemann 
problem to this location. The reconstruction step should in principle be done in characteristic variables, however for the 
second and third order schemes this can be done in primitive variables, w, as well. In the case of MHD, the latter are density 
p, pressure p, velocity v and magnetic field B. An approximation of Wij-i of an i-paritcle state at Xij is given by a first-order 
Taylor expansion from the point x^: 



Wi + ri(xjj - Xi)"(D°w)i, 



(13) 



where {D°'w)i is a gradient estimate of the primitive variables computed with Eq. [6] and is a vector of limiting functions 
which is required in order to assure non-oscillatory reconstruction ( Balsara|2004[ ) 



Ti — mm 



1, K min 



i,ngb 
^max 
i,mid 



Wi 



i,ngb 
I, mid 



(14) 



Here, w^^gb a^nd w^^g^, are the minimal and maximal w respectively over all neighbours that particles i interacts with, and 
w™^ij and w™^^;j are the minimal and maximal w resulted from the reconstruction in Eq. [13] for each of these neighbours. 
The scalar constant vector k should have values between 0.5 and 1.0 in order to achieve second order of accuracy. Following 
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suggestions of Balsara (20041, value 0.5 should be used for both pressure and velocity, and 1.0 for both density and magnetic 
field; however, we find no problems while using 1.0 for all fluid quantities. The frame velocity at Xij is approximated as 
SLij = {sLi + SLj)/2. Finally, the reconstructed left and right states, the frame velocity a^j, and the unit vector Uij are used to 
obtain the flux from the ID Riemann problem. 

Instead of a linear, a piecewise parabolic reconstruction can also be used to achieve third-order spatial accuracy. In 
Appendix|A] we describe parabolic reconstruction of a scalar field q(x, f). Due to large operation count, this reconstruction is 
presented only for purpose of completeness, and in the test that will follow later, only linear reconstruction is used. 



3 IMPLEMENTATION 
3.1 Smoothing length 

The smoothing length in our scheme is a property of the particle distribution and, in contrast to SPH, does not depend on the 
fluid state. In principle, a constant smoothing length can be used throughout the whole space and time domain. In practice 
however this cause difficulties due to the possible development of wide range in particle number densities as the simulation 
progresses. Similarly to SPH, this can lead to under- or oversampling in low and high particle density regions respectively. 
The approach used here is inspired by conservative SPH formulations (Monaghan 2002: Spring el fc Hernquist|2002 1. The idea 



is to constrain the smoothing length of a particle i, hi — /i(xi), to the particle number density at this location, rii = n(xi 



Cn.hf = iV,gb. (15) 

This tends to maintain approximately A^ngb number of neighbours for each particle; here C — 1, n and 47r/3 for D = 1, 2 
and 3 dimensions respectively, and n(xi) = l/ui(xi), where w(xi) is defined in Eq. [3] As in conservative SPH equations, hi is 
obtained by iteratively solving Eq. |15| for example via Newton-Raphson method (e.g. [Press et al.||1992 L One might be also 
tempted to use the continuity equation 

d/i(x) _ fe(x)V ■ a 
dt D ' ^ ' 

to compute time evolution of the smoothing length from its initial value. This, however, is undesirable for two main reasons: a) 
the result depends on the functional form of the divergence operator, and b) in discontinuous flows the V • a may be undefined 
at some points, which can result in unexpected behaviour. As a result, in our tests we chose to iteratively solve Eq. |15| but 
we use the differential form to predict h(x.) as a first guess to an iterative solver. 

Finally, knowledge of smoothing length permits calculation of the rest of geometric quantities, such as effective volume of 
a particle, Vi. It is possible to use numerical quadrature to evaluate J tpi{'x.) cbi with a desired accuracy, however we find that 
defining Vi — ui(xi) works fine for our purpose, and therefore we decided not to perform more accurate volume estimates. 



3.2 Particle regularity 

Particle regularity is an important aspect of the scheme. If particles are randomly sampled within a domain, there is non-zero 
probability that particle's smoothing length, h, will differ significantly from the average h in its neighbourhood. Furthermore, 
the resulting /i-distribution will not be a smooth function of position, and therefore will not be differentiable. This will 
break the approximation which lead to Eq. |1Q| Namely, the variation of h within the neighbour sphere will be large enough 
that the estimate in Eq. [4] will result in intolerable errors which produces unexpected behaviour, such as negative values 
of density or pressure. To avoid these situations, the particle distribution must be first regularised. If the initial particle 
distribution is regular, it will maintain its regularity during the simulation except in the regions where particle velocity field, 
a, is discontinuous, e.g. across shock waves ( Vila|1999 Lanson fc Vila|2008a l. The criteria which determines regularity of the 



particle distribution depends on the approximations of Eq. |4| Expanding /(x) to the first order, gives 

J /(x)i/.,(x) dx = /, J i/'«(x) dx + (V/)> ■ J{^- xOV'>(x) dx. (17) 
The first term is fiVi, and we can rewrite the integral in the second term in the following form 

(x - xOV'^x) dx = / yit;(x, + y)W{y) dy, (18) 



where in the right hand side we changed variables from x to y = x — x^. If we require that w(xi+yj)dyj ~ d is approximately 
constant in the neighbourhood of an i-particle, we can discretise the integral on the right hand side to obtain d YjV^iyj)- 
Hence, if the particles are relaxed such that this sum is minimised, the start up noise becomes negligible. Empirically, we 
found that the particle distribution is regularised if the following quantity is minimised 



^7^ = ^|AR,|^ (19) 
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where 

AR, = ^(x, -xOW(xj -x,,/iO, (20) 

Ideally, 5TZ should be equal to zero, but this is appears to be only possible if particles are arranged on a lattice, for example 
a cubic or hexagonal close-packed lattice. If particles are sampled randomly, which is more desirable in many problems, their 
positions must be adjusted until Eq. [19] reaches its (approximate) minimum before assigning fluid state to the particles; this 
will reduce the start-up noise in a simulation. Afterwards, this particle distribution can be used to assign initial conditions 
for a problem at hand. To regularise, or to relax, particle distribution, we use the following iterative procedure. First, we 
compute Eq. [20]for all particles. Afterwards, i-particle position is updated: r""*"^ — v" aAR^, where a < 0.1; such update 
reduces 5TZ. This operation is repeated until 5TZ is reached its minimum, or a desired minimal value. 



3.3 Time marching 



Due to Godunov's nature of the particle conservation laws, it is tempting to employ Hancock scheme (e.g. van Leer|2006[ ) to 
achieve a single-stage second-order accurate time integration. However, such scheme does not include a predictor for transverse 
waves in multi-dimensional Riemann problem, and therefore these remain only first-order accurate. In multiple dimensions, 
unsplit numerical schemes usually adopt Corner- Transport- Upwind method (CTU, Colella 1990 Stone et al.|2008 1 for one-step 
second-order time integration, but applicability of CTU to meshless schemes is not clear. Nevertheless, higher than the first 
order accurate time integration can be achieved with multi-stage total-variation diminishing (TVD) Runge-Kutta methods 
( Gottlieb fc Shu|1998| ). Here, we use a two-stage second order TVD Runge-Kutta time marching scheme 



(Vu)} = 1 [{VuY: + (\/u)° + {Vuf^At 



(21) 
(22) 



where (V^u)^ and (V^u)^ are time derivatives of an i-particle computed from (V^u)^ and (l^u)^ respectively via Eq. |12| Vi 
is the effective volume of i-particle, and At = cfl x min(I/i/ci,sig), where Li = Vi, ^fWipK and (3Vi/47r)^''^ for ID, 2D 
and 3D respectively is a measure of particle linear size, Csig is particle's signal speed (speed of sound for HD and of the fast 
magnetosonic wave for MHD), and cf 1 < 1 is a usual Courant-Fridrisch-Levy number. 
Particle positions obey the following equation of motion 

dxi _ 

where a; is particle's velocity. In practice we set it equal to the fiuid velocity, a^ = 
out in drift-kick-drift approach. Namely, the particles are first drifted from their current positions, x", to the position at half 
time-step 



(23) 

and the time integration is carried 



= Xi -I- ^a°Ai. 



(24) 



This particle distribution is used to compute smoothing lengths via Eq. [15] and other geometric quantities. Afterwards, we 
apply a two stage TVD Runge-Kutta method to perform an update from (Vu); to (Vu)! while keeping particles fixed in 
space and setting Vi = Vf — V^ , where V}" is i-particle volume at half time-step. Finally, the particles are drifted for another 
half time-step with the updated fluid velocity, a^ = , 



h I 1 1 A -I 

t,- H — a,- At. 
' 2 



(25) 



3.4 Non-conservative formulation 

In the case of an HD or MHD system, the conservative formulation updates total energy instead of thermal; thermal energy 
is obtained by subtracting magnetic and kinetic energies from the total energy. When the supersonic advection is present, the 
sum of thermal and magnetic energies, U — Eth + Smag, is obtained by subtracting two large numbers, namely total energy 
and kinetic energy. To avoid this, we suggest an alternative non-conservative formulation, which evolves U instead of the total 
energy. Writing E = U + P^/2M, where P and M are momentum and mass respectively, gives 

j2 \ 

(26) 



~dt ^ ~dt ~ dt \2KI 



The latter term can be rewritten as v ■ P — v M/2 resulting in the following equations for U 

dU _dE _ dP Y^dM 

~dt ~ ~dt ' ~dt ^ Y~df' ^ ' 

Here, dP/dt and dM/dt are time derivatives computed from conservative meshless equations. In this form, the total energy 
will not be conserved to the machine accuracy, but rather to the truncation error of the time-marching scheme. In other 
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words, the integration error is now lost from the system instead of appearing in the thermal energy. Furthermore, the total 
energy can now be used as a quality control indicator of a simulation. 



3.5 Modification of existing SPH codes 



The existing SPH codes which use conservative SPH formulation ( Monaghan 2002 Springel & Hernquist 2002 1 can be 
straightforwardly modified to implement our weighted particle scheme due to similarity of SPH equations of motions and our 
equations of meshless conservation laws, Eq. |12[ In particular, the neighbour search should be modified such that Eq. |15| is 
solved, which constrains the number of particles, instead of the enclosed mass, in the neighbour sphere. At the end of this 
step, the hi for each particle will be known that permits to compute volume of a particle, Vi = Wiihi). This volume is required 
to convert conservative fluid variables, (V^u)i = ViUi, to primitive ones Wi. Afterwards, the first loop is carried out over 
gather neighbours which computes the renormalisation matrix, Eq. [T] and gradients of primitive fluid variables. The second 
neighbour loop is required to compute the limiting functions, Eq. |14[ and this loop must be carried out over both gather and 
scatter neighbours due to need to limit reconstruction to each of the interacting particles. Finally, in the third neighbour loop 
the interactions between the particles are computed. This is done in exactly the same way as in SPH, except that for every 
j-neighbour of an z-particle, the fluid states are reconstructed at the midpoint, x^j = (xi + Xj)/2 through Eq. |13| Finally, 
these together with the vector n^j in Eq. |12| and the interface velocity a^j = (v^ + Vj)/2, are used in the Riemann solver to 
compute the interface fluxes, (Fij • iiij), which we describe in the following sections. At the end of this final neighbour loop, 
one will have time derivative that should be used to update the conservative fluid variables. The global conservation of mass 
and other conservative quantities, in the absence of source terms, is maintained to the machine precision independently of the 
particle distribution. This also includes total energy, unless Eq. [27] is used, in which case the total energy is conserved to the 
truncation error of a time integration scheme. 



4 FLUID DYNAMICS 

Our meshless conservative equations can be applied to a system which can be written in the form of conservation laws. 
The physics in this case is completely described by the source terms, S, and fluxes, J-. The latter can be obtained by 
solving an appropriate one-dimensional Riemann problem between a particle and its neighbours. Hence, the Eq. [12] can 
be applied to a variety of problems, such as hydrodynamics, magnetohydrodynamics, and radiative transfer in the flux- 
limited diffusion approximation. In what follows, we present application of our scheme to both ideal hydrodynamics and 
magnetohydrodynamics. 



4.1 Ideal hydrodynamics 

The Euler equations of ideal hydrodynamics in a frame moving with the velocity a read 

dU 



dt 



+ V- (J^-aW) = 5, 











etot 













(28) 



where 



(29) 



where I is a unit tensor, and other symbols have their usual meaning. 

To obtain fluxes, an ID Riemann problem is solved between two particles in the n^j direction. This is accomplished by 
deflning the rotation matrix, such that in the new coordinate system n^j coincides with the x'-axis, i.e. n^j = An.ij = 
(|nij| , 0, 0). This transformation is applied to all vector quantities from both the left and the right states, the scalar quantities 
are left untouched. For example, the velocity transformation results in = A^k = (t'i.xi '^j/./f ) ^z,if )i where K = L for the 
left and K = R for the right state. These transformed states are the initial conditions of the Riemann problem for ID Euler 
equations 



dt 



dx 



0, 



(30) 



where, G'lo = J^'id — a'^U'io and 



P 

etot 
pv'x 

pV'y 

pv'z 



\ P< } 



( P-V'x \ 

(etot + p)v'x 
pvWx +P 

pvWx 



(31) 
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Figure 1. This figure shows the wave-structure of HLLC Riemann solver. Namely, HLLC solver resolves all three HD waves: left, right, 
and middle wave also known as contact discontinuity waves. The left and right wave can be either shock or rarefaction wave, or both. 
This solver approximates the structure of the two intermediate states, between Sl and Sm, and Sm and Sn, by a constant state. 



The space discretisation of this equation has the following form 



dt 



+ 



Q',. 



Ax 



0, 



(32) 



where Ax is size of the grid cell in ID, and ^/i4-i/2 is an interface flux between cells i and i + 1. This is the flux required to 
substitute into Eq. |12[ and it can be obtained from the solution of an appropriate ID Riemann problem. 

The solution of this ID Riemann problem gives the fluxes, Q' , for each of the component of U' in direction Uij. While the 
fluxes of the scalar components can be directly used in Eq. |12| those of a spatial vector, however, must be rotated back to 
the original coordinate system because the flux of an a;-component of the vector, rather than x' , is required in the direction 
iiij. In the case of the velocity vectory, one first combines (f^, v'y,v'^)-iiux into a vector G'(v') = {Q'{v'^), Q'{v'y), G'iv'^)). Then 
an inverse transformation is applied to compute {vx,Vy,Vz)-&ux, namely G'(v) = ^~^G'(v') = (Q' {vx),G' ii'y),g' (vz)), and 
these fluxes are then substituted into Eq. |12| 

This approach permits the use of Riemann solvers that directly approximate interface flux, rather than fluid states. In 
what follows, the HLLC Riemann solver is used to compute the flux. The HLLC solver requires only velocity estimates for 
certain characteristic waves, and is oblivious to the exact form of the equation of state. As a result it can be easily extended 
to equations of states other than that of ideal gas. Here, we present formulae of the HLLC Riemann solver in moving frame 
which can be di rectly used in Eq. |12| for the detailed derivation we refer the interested reader to published literature (e.g. 
Toro|1999| (§10), |Miyoshi fc Kusano|2005l ). 

A typical wave-structure of HLLC solver is shown in Fig. [l] The velocity of left, right and middle waves are Sl, Sr and 
Sm respectively, and the constant states sandwiched between these waves are U'^ and Ur respectively. The HLLC-flux at the 
interface moving with velocity a'^ is 



T'r + Sr{Ur - U'r) - a'xUR 



a'x < Sl, 
Sl < a'x < Sm, 
Sm < a'x ^ Sr, 
Sr < a'x. 



(33) 



where — T{U'l) and F'r = F{U'r). The intermediate states are defined by 



Sk - v'xK 

Sk — Sm 



Ek 
pk 



{Sk 



1 

Sm 

v'zK 

:) {Sm 



+ 



Pk 



pk(Sk--"^k) J 



(34) 



for K = L and K 
Cs = max(csL,CsH) 



— R. Finally, estimates of wave speeds are Sl ~ mm{v'x l,v'xr) — Cs and Sr = max(w^i, w^jj) -I- Cs, where 
Other wave-speed estimates, such as those based on Roe-averages, can be used as well. Finally, the speed 



of the middle wave and the pressure in the ★-states is given by the following formulae 

r. _ PR-PL+ Pv'xl{Sl - v'xl) - Pv'xr{Sr - v'xR 



Pl{Sl 



Pr{Sr 



(35) 



p = 



{Sr — Vxr)prPl — {Sl — Vxl)plPr + PlPr{Sr — Vxr)(Sl — Vxl){vxR — Vxl) 
{Sr — ur)pr — {Sl — ul)pl 



(36) 



with vt 



- VxR = Sm and pi = p*r 
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4.2 Ideal MHD 



Our meshless conservative equations can also be applied to the equations of ideal MHD. In particular, we can rewrite MHD 
equations in the conservative form with the following conservative variables, fluxes and source terms 



P 
etot 
pv 
B 



(etot + Pt)v ~ (v ■ B)B 
pv ® V + PtI - B (g) B 
V ® B - B (g) V 



(37) 



V / 



here Pt = p + B'^ /2 is sum of the thermal and magnetic pressures. The biggest difficulty in solving these equations is to 
maintain V • B = constraint. While there are several ways to satisfy this constraint in finite-difference methods, it is not 



clear how this can be done in meshless schemes. In SPH, Rosswog & Price (20071 were able to circumvent this problem via 



use of the Euler potentials, a and /3, which in smooth flows are advected with the flow, and are used to compute magnetic 
field B = Va x V/3. While mathematically this guarantees zero-divergence, in practice, however, the divergence is non-zero 
because of non-commuting nature of SPH cross-derivatives. This is also holds for the vector potential. A, which may explain 
unstable behaviour of SPMHD equations with vector potential ( Price|2010 1. While SPH equations appear to be stable in the 
Euler potential formulation, they are unable to model topologically non-trivial field configurations ( Brandenburg|[20T0 1 , and 
therefore are not suited to simulate complex magnetic field evolution, such as winding of magnetic field lines and magneto- 
rotational instability. Motivated by this, and by the fact that the existing MHD Riemann solvers use magnetic field as a 
primary quantity, we chose to evolve B instead, and apply one of the divergence cleaning methods known to work in Godunov 
finite-difference MHD schemes. 

In their paper, [Powell et al.| ( |1999| ) showed that a self-consistent MHD equations must include source terms proportional 
to V ■ B to insure stability and Galilean invariance 



-V B 



V B 
B 

V V / 



(38) 



While this formulation both removes the instabilities associated with non-zero divergence and tends to keep the divergence at 
the truncation level, the divergence can still grow in certain situations. To avoid this, we also include a Galilean invariant form 
of the hyperbolic-parabolic divergence cleaning method due to |Dedner et al.| ( |2002[ ), which results in the following system 



f P \ 



U 



P 
etot 
pv 
B 



T ■ 



(etot + Pt)v - (v • B)B 
pv ® V + PtI - B ® B 
V (g) B - B (g) v 
V pV'v / 



/ 

-(V-B)v 
-(V-B 
-(V-B 
V -(V-B 



B B Vt/) 



B 



ChP 



- 7/;p/r 



(39) 



Here, Ch and r is the speed and the damping timescale of the divergence wave, also known as 8-wave, respectively. Usually, 
Ch is set to be highest characteristic speed, i.e. that of the fast magnetosonic wave, r = L/^CtCh), where L is an effective size 
of the particle defined in i3.3 and Cr is a constant, which, following Mignone & Tzeferacos (20101, we set Cr 



0.03. 



In their paper, Dedner et al. ( 2002 1 also presented a method to generalise an arbitrary ID MHD Riemann solver to 



include an additional scalar field ^j) which transports the divergence away form the source and damps it. In general, the left 
and right reconstructed states in coordinate system A (see S 4.1 1 have discontinuous B^-the magnetic field component normal 
to the interface. The ID Riemann solvers, however, require this field to be continuous, and this can be computed with the 
following equations 



B'x = 7,iB'xL + B'^r) 



2 



; {ipR - iPl) 

-h.ij 

{B'xr — -B^l)- 



I Ch, 



(40) 
(41) 
These 



Here, B'^f^ are the left {K = L) and right (K — R) reconstructed states for B'^ and tp, and Ch,ij ~ max(ch_; 
interface values of B'^ and are used to compute V • B and Vtp respectively, which are required for the source terms in Eq. |39[ 



y,(V-B). = -^Bi|n,,| 

j 



(42) 
(43) 



In Appendix [B] we provide standard formulae to compute HLL and HLLD fiuxes for the ID MHD Riemann problem, which 
take B'j. as the continuous normal component of the magnetic field. Finally, the advection flux for the scalar p^p is 



Fip Fp X ■l^ otherwise, ^^^'^ 

where, Fp is a mass flux given by the Riemann solver. This completes the description of our meshless MHD scheme. 
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5 SCHEME VALIDATION 



The weighted particle scheme is validated on several standard problems for ideal hydrodynamics and MHD problems, which 
are usually used to test hydrodynamic and MHD schemes (e.g. T6th|[2000 Stone et al.|2008 In all simulations what follows, 
we use TVngb = 19 and 32 in Eq. [15] for 2D and 3D simulations respectively. This choice is motivated by the analogy with 
finite-difference schemes. In one dimension, the number of the neighbouring cells that a given cell interacts with depends on 
the order of the numerical scheme, and is usually two for the second order scheme. In SPH, two to four neighbouring particles 
are usually used in one-dimensional simulations. If this number is scaled to three dimensions, and taking into account that 
a kernel has spherical shape, the estimated number of neighbours are 16 and 32 for 2D and 3D respectively. The MHD 
simulations in 2D appear noisy with A^ngb = 16, which motivated us to use larger A^ngb. In principle, large A^ngb can be used, 
should this be necessary, but this will decrease the resolution due to larger smoothing length. 



5.1 Shock tubes 

5.1.1 Brio-Wu shock tube 



This problem was introduced by Brio & Wu ( 1988 1 to test the ability of an MHD scheme to accurately model shock waves, 
contact discontinuities and compound structures of MHD. Here, the problem is solved in a periodic 2D domain of size 
[0, 4] X [0, 0.25] with randomly sampled 5 • 10* particles; this results in an effective resolution of 895 x 56. The particles are 
initially relaxed before the initial conditions are set. The left state, a; < 2, is set with the following values: — 1, Pl ~ 1, 
ByL ~ 1. The right state has pu — 0.125, pr — 0.1 and ByL = —1. Both states have zero initial velocities, Bz = and 
Bx = 0.75. The problem is solved with an ideal gas equation of states and 7 = 2. Various profiles, e.g. density and pressure, 
at time t = 0.2 are show in Fig. [2] 

The results are overall consistent with ID Godunov-MHD scheme, namely jumps across the discontinuities and location 
of discontinuities. The two bottom right panels show the value of B^ field and the divB — LV ■ B/jB| as a function of particle 
x-coordinate. The parallel magnetic field slightly deviates from its constant value, except near discontinuities where it exhibits 
jumps. The divergence, however, remains small, even across discontinuities. The existence of blip in pressure at location of 
contact discontinuity, x ~ 2.1, and shock waves, x ~ 2.3 has the same origin as in SPH: the particle distribution across a 
discontinuity is less regular in a sense that the approximation in Eq. (|4]is not sufficient to provide accurate results. This is a 
known issue in Godunov SPH, and higher order approximations are able to reduce the amplitude of the blip ( |Inutsuka|2002[ ) . 
Overall, the solution obtained by the meshless scheme is in a good agreement with high-resolution ID Euleriean schemes. 
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5.1.2 Toth shock tube 



Another challenging shock tube problem was introduced by Toth (20001. 
supersonically collide with each other. 

V • B = produce wrong jump conditions. This was challenged b y [Mignone fc Tzeferacos] ( |2010[ ), who showed that if 



In this problem, two streams of magnetised gas 
Toth (2000) showed that some MHD schemes with source terms proportional to 



divergence is cleaned in hyperbolic-parabolic manner (Dedner et al. 20021, the jump conditions are correct even with the 



presence of the source terms. The problem is set in 2D periodic domain with size [0, 4] x [0, 0.25] with randomly distributed 
5 • 10* particles. After the particle distribution is relaxed, the following initial conditions are set. The left states has pL ~ 1, 
PL = 20, VxL ~ 10, and the right state has pn = 1.0, pr = 1, v^r = —10. Both states have Bx = By = 5/\/47r. This 
problem is solved with an ideal gas equation of state and 7 = 5/3. The solution to this problem is shown in Fig.|3] which 
can be compared to solutions obtained by (| T6th|[2"000[ ) and 'Mignone fc Tzeferacos| ( |2010[ ). Our meshless scheme is able to 
recover correct jump conditions and to maintain constant Bx field within few percent accuracy, except across discontinuities. 
Furthermore, the divergence in this problem remains small. 



5.2 Advection of a magnetic field loop 

This problem tests the ability of the scheme to transport magnetic loop across computational domain. This problem is 
proven to be a stringent test for finite-difference schemes. The computational domain in this test is a cuboid with dimensions 
[0, 2] X [0, 1] X [0, 0.5] which is initially filled with 128 x 64 x 32 particles on cubic lattice to assure zero noise. The fluid 
state is set to p = 2 inside the loop and p = 1 outside the loop, p = 1, -v — (2,1,0.5) and B = {f(R)y,—f{R)x,0), where 
R = ^ x'^ + y'^ and 

f(F!] — l -^0/^ < R < Rf), , > 

^ ' \ 0, otherwise. ^ ' 

Here, Ro = 0.3 is a radius of the loop. Bo — 10~^ is the initial magnetic field strength that results in 2/3 = Pgas/^mag ~ W^. 
With such high /3 magnetic field does not play dynamical role and should be transported as a passive scalar. Periodic boundary 
conditions are used in this test, and the ideal gas equation of state with 7 = 5/3. 

In Fig. [4]we shows magnetic field structure at the t = and t = 10, which corresponds to ten crossings of the computational 
domain, and in Fig. [5]we show the magnetic energy as a function of time. This figures demonstrate the ability of the meshless 
scheme to advect magnetic loop quite well. Furthermore, the decay of the magnetic energy is at least as slow as resulted 
from high-order Godunov MHD schemes. This is certainly expected in light of semi-Lagrangian nature of the scheme. One 
may however expect that the energy should not decay at all since the magnetic field is transported as a passive scalar. 
Indeed, Fig. [6] demonstrates that the scheme is able to transport mass, and therefore passive scalar, without any diffusion. 
The difference with magnetic field stems from the different nature of equation that transport mass, or advect passive scalar, 
and the induction equation as implemented in the scheme. In fact, the decay is caused solely by the diffusion of magnetic 
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Figure 4. Amplitude (colours) and direction (vectors) of magnetic field at the beginning (left panel, t = 0) and at the end (right panel, 
t = 10) of the simulation. The magnetic field distribution is show in the XY-plane passing through z = 0.25. In the right panel, the 
strength of magnetic field in the centre is close to zero due to numerical resistivity, which is consistent with finite-difference calculations. 



LU 




Figure 5. This figure demonstrates decay of magnetic energy during advection of magnetic field loop. The magnetic energy decay is due 
to numerical resistivity introduced by divergence cleaning procedure, and not through the solution of an MHD Riemann problem. 



field due to divergence cleaning equations, i.e. Eq. (40 and Eq. [41] The fiuxes resulting from HLLD Riemann solver are in 
fact zero since both particles and the frame move with the same velocity. The divergence of the magnetic field, even though is 
zero analytically, is not necessary zero in the discretisation set by Eq. [40] andEq. |42| and this produces evolution of magnetic 
field due to non zero value of VV' in the source terms of Eq. [39] Among all discretisations studied by Toth (20001, this is 
the special one because of its use in the discretisation of Maxwell stress term to obtain Lorentz force. If the V • B = in 
this discretisation, no force parallel to magnetic field exists, and therefore the source terms proportional to this divergence 
vanishes. Incidentally, this is the discretisation of divergence that is enforced to zero in constrained transport (CT) formalism 
( Evans fc Hawley][1988 K More importantly, the maintenance of zero divergence in other discretisations, such as cell-centred, 
does not guarantee vanishing divergence in CT-discretisation, but this will probably be small thus giving minimal damage to 
the solution. 
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0.2 



Figure 6. This figure sliows tliat our scheme is able to transport mass, and therefore passive scalars, without any numerical diffusion 
in a constant velocity field. The right and left panels show plots of the density field in the XK-plan passing through z = 0.25, at the 
beginning and the end of the simulation respectively. 



5.3 Blob test 



An interesting and challenging problem to test particle hydrodynamic schemes has been proposed by Agertz et al. (20071. This 



problem demonstrates the destructive property of the Kelvin-Helmholtz instability (KHI) in three-dimensional simulations. 
The setup consist of a dense cloud with pci = 10 moving supersonically through a less dense ambient fluid with Pamb ~ 1- 
Initially, the cloud is at rest and in the pressure equilibrium with the ambient fluid, at Pamb ~ 1. The velocity of the ambient 
medium is v = (2.7cs,amb, 0, 0), where Cs,amb is its sound speed. The cloud radius is Rd — 0.1, and an ideal gas equation 
of state is used with 7 = 5/3. The initial magnetic field is set to zero, and it will stay so throughout the simulation. This 
problem is set in a periodic domain with dimensions [0,3] x [0, 1] x [0, 1], in which 7 • 10"" particles are sampled in the strip 
\y — 0.5| < l.li?ci and 3 ■ 10^ outside. This setup permits to resolve cloud and the impacting ambient fluid with high reso lution 
{h ~ 0.01) for a total of 10^ particles. The particles are sampled from a three-dimensional Sobol quasi-random sequence (Press 



et al.|[l992 l. To remove the initial noise, this initial particle distribution is relaxed before the initial values for the density 
field, pressure and velocity are set . 

Following Agertz et al. (20071, the Kelvin-Helmholtz timescale is defined Tkh ~ 1.6rcr, where r^r — 2Rc\ \/ pd /Pamb / «amb 



For parameters used in this simulations, this gives Tkh ~ 0.3. The snapshots of the density distribution in plane z = 0.5 at 
t = 0.5, 1.0, 1.5 and 2.5 Tkh are shown in Fig. [7| These are in an excellent agreement with those presented in [Agertz et al.| 
( |2007[ ). The time dependence of cloud mass is shown in Fig. [8] where a particle is considered to be part of a cloud if p > 0.64pci 
and T < O.OTamb. In agreement with finite-difference methods, the cloud lost nearly 90% of its mass within 2 Tkh. 



5.4 Spherical blast-wave 

The problem is initiated with an overpressured central region in a uniform density and magnetic field. The computational 
domain is a unit square filled with fluid with p = 1. Within R < Ro = 0.1, the pressure is set to 10, whereas outside p = 0.1. 
In magnetised case, there is also a uniform magnetic field B — (l/\/2, l/\/2,0). The equation of state is that of an ideal gas 
with 7 = 5/3. The particles are sampled in a periodic box [0, 1] x [0, 1.5], such that 5 ■ 10* particles are randomly sampled 
in three nested rectangles: [0,1] x [0.1,5], [0.25,0.75] x [0.375,1.125] and [0.375,0.625] x [0.5625,0.9375]. Before the initial 
conditions are set, the particles are relaxed into a regular distribution to reduce start-up noise. In figures |9] and [lO] we show 
density plots for non-magnetised and magnetised cases respectively. Of particular interest here is the ability of the particle 
weighted method to resolve Richtmyer-Meshkov instability, shown in right panel of Fig. [9] In the magnetised case, however, 
the presence of strong magnetic field inhibits development of this instability (the right panel of Fig. |10[). 



5.5 Orszag-Tang vortex 



The Orszag-Tang vortex ( Orszag & Tang 1979 \ is a standard test problem that is used to validate many numerical MHD 
schemes. The setup involves periodic domain of size [0, 1] x [0, 1] with an adiabatic equation of state with 7 = 5/3. The 
initial density and pressure are set in all computational domain to 25/(367r) and 5/(127r) respectively. The velocity v — 
(— sin(27r2/), -I- sin(27ra::), 0) and magnetic field B — (— Bo sin(27rj/), _Bo sin(47ra;), 0.0), where Bo = l/x/ivr. The second simulation 
involves the same initial conditions, with the exception that the problem is solved a boosted frame, with the initial velocity 
V = Vrcst + Vboost, where Vboost = (10, 10, 10). The simulation is solved with 10^ particles, where first 5 • 10* where randomly 
sampled within the square [0, 1] x [0, 1] and the second 5- 10* in the square of half size, [0.25, 0.75] x [0.25, 0.75]. This permitted 
to achieve high resolution in the central region of the problem. The initial conditions were set, as soon as distribution was 
relaxed. 
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Figure 7. Density field plot in XY-plane passing through z = 0.5. The panels show density profile at time 0.5, 1.0, 1.5 and 2.5 Tkhi 
where Tkh ~ 0-3. At t = 0.5Tkh (top panel), the cloud shape is distorted by RAM pressure. At t = Tkh (second panel from the top), 
the KHI instability deforms the shape of the cloud further. The cloud destruction begins at i = 1.5Tkh, and hy t = 2.5Tkh the cloud 
lost more than 90% of its mass. 

In the Fig. [Tl]we show density at t = 0.5 for both rest-frame and boost-frame initial conditions, and in Fig. |12| we 
show particle distribution. As expected, there are no discernible differences, and both simulations resolve discontinuities well. 
Furthermore, due to Lagrangian nature of the method, the particle number density correlates with the mass density. It is 
worth pointing out that in this particular simulation, the particle distribution appears to remain regular even in the vicinity of 
the shocks. The Fig. |13|show ID pressure profile at y = 0.3125 (top) and y — 0.427 (bottom) for both boosted and rest-frame 
simulations. The agreement with finite-difference scheme is excellent, ac can be compared to published results (e.g. |T6th|2000 
[Rosswog &: Price|2007| [stone et al.|2008| |. In contrast to shock-tube problems, no pressure blips are visible here. 



5.6 MHD rotor 



The rotor problem, introduced by Balsara & Spicer ( 1999 1 to test propagation of strong torsional Alfven waves, is also 
considered as one of the standard candles to validate numerical MHD schemes. Here, the computational domain is a unit 
square, [0, 1] x [0, 1]. The initial pressure and magnetic field are uniform with values p = 1 and B = (5/\/47r, 0, 0). Inside 
R < Ro ~ 0.1 there is a dense uniformly rotating disk with p = 10 and v = (~2y / Ro, +2x/ Ro,0), where R = + y^- 
The ring Ro < R < Ri = 0.115 is occupied by a transition region with p^l + 9f{R) and v = {-2yf{R)/R, +2xf{R)/R,0), 
where 



fiR) = 



1 

Rl-R 
Hi -Ho 





R < Ro, 
Ro<R<Ri 
Ri < R. 



(46) 



Outside R > Ri, the velocity is set to zero and p = 1. This problem uses an ideal gas equation of state with 7 — 1.4. The 



particles are distributed in the same way as in strong blast wave problem (S5.4I. The problem is solved in both rest and 
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time [Tkh] 

Figure 8. This plot sliow cloud mass as a function of time. In agreement with finite-difference scheme the cloud looses most of its mass 
at t ft! 2 Tkh I and is completely destroyed at t > 3Tkh- 




Figure 9. This figure show the density at t = 0.2 and t = 1.5 for non-magnetised spherical blast wave problem. The location of the 
shock front is in agreement with high-resolution conservative numerical schemes jStone et al.|2008l l. Right panel show dense fingers in 
rarefied media which are formed by Richtmyer-Meshkov instability. This demonstrates that our scheme is able to capture important fiuid 
instabilities without any fine-tuning. 
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Figure 10. This figure show the density at t = 0.2 and t = 1 for magnetised spherical blast wave problem. The initial magnetic field in 
this problem has angle of 7r/4 with x-axis. Since the shock moves more easily along the magnetic field lines, this explains the elongated 
shape of the shock-front, which contrasts with hydrodynamical case. As in hydrodynamical test, here the location of shock from is 
in excellent agreements with conservative Eulerian MHD schemes | |Stone e"t~ al. 2008). The right panel, shows further evolution of the 
black wave, after shock reached the rarefied medium. In contrast to the hydrodynamical case, magnetic field inhibits development of 
Richtmyer-Meshkov instability. 




Figure 11. Density distribution of Orszag-Tang vortex at t = 0.5. The left and right panels show the result of the simulation in the rest 
and boosted frame respectively. We point out that even at low resolution our scheme is able to resolve the fine structure in the bottom 
left quadrant; in this quadrant, the effective resolution is 128 X 128, compared to 256 X 256 in the centre. 
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Figure 12. Particle distribution of Orszag-Tang vortex at t = 0.5. Tlie left and riglit panels show the particle distribution in all domain 
and in the central region respectively. Due to Lagrangian nature of the simulation, the particle distribution correlated with that of the 
density. We would like to note that even across the shocks, the particle distribution remains remarkably regular without any observable 
signs of clumping. 




Figure 13. Pressure distribution of Orszag-Tang vortex a.tt = 0.5 and y = 0.3125 (top) and y = 0.427 (bottom). This can be compare to 
the results of other finite-difference or meshless MHD schemes (e.g. [Rosswog fc Price|2007[ [Stone et al.|2008[ |. In particular, our scheme 
has little noise in pressure, and the shock waves are resolve within few particles, without observable oscillatory behaviour in post-shock 
regions. 



boosted frames, with iiboost ~ (5/0.15,5/0.15,5/0.15). The plots of density, gas and magnetic pressure, and mach number 
are shown in Fig. |14[ The ID slices of magnetic field at a; = 0.5 and y = 0.5 are show in Fig. ^] to facilitate comparison to 
Eulerian MHD schemes. 



5.7 Magneto-rotational instability 

Magneto-rotational instability, or MRI for short, is a powerful local instability in weakly magnetised disks that plays an 
important role in astrophysics ( Balbus fc Hawleyl|1991 19981. The ability of the scheme to model this instability is crucial 
for simulation of magnetised accretion disks, or any other simulation where non-trivial MHD efi'ects are expected to appear. 
Below we conduct two simulations that test the ability of the meshless MHD scheme to faithfully model MRI. 
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Figure 14. Density (top left), gas pressure (bottom left), magnetic pressure (top right) and Mach number (bottom right) plots for MHD 
rotor problem at t = 0.15. The Mach number panel plots Mach number using rest-frame velocities, that is first subtracting the boost 
velocity and afterwards computing mach number. The result of the rest frame simulation have little difference, and therefore not presented 
here. The colour coding is chosen to facilitate the comparison with published results produced by FLASH3 code for this problem | |Fryxell| 
|et al.|2000| Sect 21.2.3 in FLASH 3 user guide manual ,http : //flash ■ uchicago ■ edu/website/codesupport/f lash3_ug_beta^ . 



5.7.1 2D axisymmetric shearing box 

The first simulation is solved in a local 2D axisymmetric shearing sheet with the initial conditions exactly the same as in the 
fiducial model of |Guan fc Gammie| ( |2008| . To repeat, the unit box has a size [0, 1] x [0, 1] with an initially uniform density fluid 
with unit density, po = 1, embedded in a vertical magnetic field Bz = Bo sin(27r2;/AMRi), where Amri = 27r-\/l6/15uA/f2 is the 
fastest growing MRI wavelength, va ~ Bq/ y/po is Alfven speed, and Q is angular velocity which is set to unity throughout the 
simulation. The magnetic field strength Bq = y^2po7A), where Po = 1348 in order to excite m = 4 MRI mode. The instability 
is seeded by random velocity perturbation Sv = O.OlCs. Isothermal gas equation of state is used, p = c^p, with Cs = 1 in all 
domain at all times. The boundary conditions are periodic in z-direction, shear-periodic a;-direction, i.e. 



f{x, z) = f{x + nj:Lx, z -I- n^L^) 



(47) 



Vy(x, z) = Vy{x + rixLx, 2 -I- rizLz) + nj:qi}Lx. 



(48) 
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Figure 15. Distribution of By in y = 0.5 slice (top pannel) and Bx in x = 0.5 slice (bottom) for MHD rotor problem at 
plots can be readily compared to Fig. 26 of ( [Stone et al.|2008| . The agreement is rather good everywhere, except for By 
and X 0.66, where the meshless scheme appears to diffuse the discontinuity in magnetic field. 



t = 0.15. This 
near x ^ 0.34 



Here, q = —l/2dlnQ/d\nR which is equal to 3/2 for a Keplerian disk, and 
shearing box model, the momentum equation have the following form 



and Uy are arbitrary integer numbers. In this 



dtipv) = - V ■ (/5V ® V + PtI - B ® B) - 2n X (pv) + 2pqVrxx, 



(49) 



where Pt = p + /2 is the sum of the thermal and the magnetic pressure, the second and third terms on the right hand side 
are Coriolis and centrifugal forces respectively. 

This problem is solved in a periodic domain, with 128 x 128 and 256 x 256 particles initially distributed on Cartesian grid. 
The toroidal, By, component of magnetic field is show in Fig. [16] at f = 5, 10, 15 and 20. The m — A MRI mode is clearly seen 
at t = 10, and the bottom right panel demonstrates the break up of the laminar flow into turbulent. Fig. [17] shows both the 
toroidal magnetic fleld (left panel) and divB (right panel) at t = 40. The aim is to demonstrate that even in turbulent flows, 
the scheme is able to keep the divergence small. Time dependence of magnetic field energy and Maxwell-stresses is show in 
Fig. |18[ The initial growth can be fit with Emug/U oc exp(0.75t), and is in excellent agreement with Fig. 4 and Fig. 11 of|Guan 



& Gammie (20081. The decay at late times differs, which is an expected result, since this depends on the details of numerical 



dissipation, which differs among MHD schemes. 



5.7.2 3D Global disk simulation 



The final test problem studies the evolution of a circular disk around a point-mass gravitational source. The aim is to test 
the particle scheme on realistic astrophysical simulation, and to determine its weak points. The problem is set up in a 



computational domai n of size [0, 20] x [0, 2 0] x 
Following methods of 
from its inner edge, R 



Lodato & Rice 



(20041 and 



0, 20] with a gravitatio nal source of unit mass located at origin (10, 10, 10). 

particle 

10)-^ is the distance from the midplane 



Alexander et al. 



( 2008 1, a circular stratified disk is sampled with 10 



1, to its outer edge, Rout = 4, where R= ^^{x — 10)^ + {y 
of the disk to the gravitational source. The scale height of the disc is H/R = 0.1, and only one scale height is sampled to avoid 
large density variations. Particles which fall within -Rmin = 0.25iiin and outside the computational domain are removed from 
the system. The particle distribution is regularised to minimise the start-up noise before the initial conditions are assigned. 
Since this problem uses open boundary conditions, we also imposed a lower and upper limit to the number of neighbours; 
specifically, we set A^ngbjow = 8 and A'ngb,up = 128 for the lower and upper limit respectively. We did this in order to prevent 
particles close to the boundaries to have either too few or too many neighbours. The lower limit, however, not been reached 
in our simulations. 

The initial density in the disk is p{R, z) — Pmid exp{—z^/2H^), where pmid = 1 is the density in the mid plane. The scale 
height is H = Csnd/'^circ^i, where Vdic = \/l/-R and Csnd is sound speed of a particle. The latter is set to be constant in time, 
but has the following spatial variation Csnd ~ {H/ R)vcirc. Pressure is set via isothermal equation of state, p = c^ndP- The 
gravitational acceleration is split into horizontal and vertical component, to make sure that the above hydrostatic equilibrium 



is satisfied, namely Og 



g{R){x - 10,y - 10, z - 10), where g(R) 



-1/R . Two simulations where run with above 
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Figure 16. Toroidal, By, magnetic field at t = 5 (top left) , t = 10 (top right), t = 15 (bottom left) and t = 20 (bottom right) for 2D 
shearing box simulation at 256 X 256 resolution. At t = 10 panel, the m = 4 MRI mode can be seen, in an excellent agreement with the 
MRI theory, and at t = 20 the flow becomes turbulent in agreement with finite-difference MHD schemes. 



initial conditions, one non-magnetised and one magnetised. In the latter case, initial constant vertical magnetic is set with 
/3 = 1348 in the midplane, which have following _R-variation, B{R) = Bosm{2n{R — 2)), and the magnetic field outside the 
ring 2 < i? < 3 is set to zero. With this setup, the orbital period of an inner disk orbit is 2n. The magnetic field is forced to 
be zero for R < 1.2 and R > 6.5. 

Density structure of the disk in XZ-plan passing through the gravitational source is shown in Fig. [19] for t = 10, 50, 100, 
150. The right panels on the figure show density distribution for non-magnetised case. As expected, the flow remains laminar 
throughout the simulation. However, at t = 100, the inner edge of the disk is noticeably damaged due to particle loss through 
Rmin- This error, generated at the inner edge of the disk, slowly propagates outward, as can be seen by the further damage 
at t = 150. This example demonstrates the importance of proper boundary conditions in meshless scheme. However, it is no 
clear how to define these. Nevertheless, the global structure of the disk remains consistent with the initial conditions and 
vertical hydrostatic equilibrium through the simulation. 

In the left panel of the Fig. ^] we show density distribution of the magnetised case. In contrast to non-magnetised 
simulation, the laminar motion breaks into turbulence at t ~ 100. For earlier times, the laminar flow permits the magnetic 
field to grow to (/3^^) ~ 0.1 as shown in the top panel of the Fig. |20[ Of the particular interest here, is the value of the 
magnetic stresses, qm = — BrBe/-Pgas, which determines the accretion rate in the disk. Time dependence of the volume 
average magnetic stress is shown in the bottom panel of the Fig. |20| The growth continues until t « 100, after which the 
flow become turbulent and the volume averaged magnetic stress remain roughly constant at (um) ~ 0.01. In contrast to 2D 
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Figure 18. Time dependence of volume averaged magnetic energy and Maxwell stresses as a+ function of time. This figure can be 
readily compared to Fig. 4 and Fig. 11 of Guan & Gammie (2008). The initial growth rate is in excellent agreement with their results. 
The subsequent turbulent evolution is different, however. This is a known result, since the turbulent evolution of the flow depends on 
the numerical scheme employed. Due to anti-dynamo theorem, the magnetic energy and Maxwell stresses decay in this 2D simulation. 



axisymmetric shearing sheet simulation, neither magnetic energy nor magnetic stressed decay with time, implying the dynamo 
activity in the disk. 

The density distribution in the mid plane of the disk is shown in Fig. |21| The left and right panels of the figure display 
non-magnetised and magnetised cases respectively, and the top and bottom panel show the profile a,t t — 100 and t = 150. The 
turbulent structure of the magnetised disk for t > 100 is apparent through the existence of small scale structures in density, 
whereas the non-magnetised disk remains laminar, with the exception of small spiral waves which are caused by the error 
propagating from the inner disk boundary. The magnetic field structure is shown in Fig. [22] for t = 50 (top panel) and t = 100 
(bottom panel). The left panel shows the amplitude of the toroidal magnetic field. Be, and the value of Maxwell stresses, 
Br Be- At t = 50, the structure begins to appear, but the magnetic field is weak to substantially influence the dynamical 
evolution of the disk. At t = 100, however, the magnetic field is amplified to a large values via MRI mechanism, and has 
non-negligible effect on dynamics. Furthermore, the toroidal magnetic field exhibits reverse, cause by the shear, in agreement 
with the theoretical expectations. 

Overall, the growth of magnetic field is consistent with published 2D and 3D shearing box simulation. Namely, the linear 
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Figure 20. This figure shows the growth of the volume average magnetic field and Maxwelian stresses in the magnetised disk simulation. 
In agreement with the 2D shearing sheet simulation, the magnetic field and stresses reach its maximum after about 15 inner orbits (or 
5.6 orbits at JJ = 2. In contrast to 2D case, however, the magnetic field does not decay furthermore, but is maintained by dynamo action. 



MRI regime is able to amplify magnetic field to (/3~^) > 0.1 until it breaks up into turbulence. Due to limited numerical 
resolution of this simulation, we were unable to resolve fastest growing MRI mode, which explains the slower than expected 
growth. Nevertheless, the magnetic stresses are roughly ten percent of magnetic energy, (om) ^ 0.01, and drive the accretion 
of the matter. This matter is accumulate at _R < 1.25 where magnetic field is set to zero at boundary condition, and explains 
dens blob of matter in the two lowest left panels of Fig. [19] 



6 DISCUSSION AND CONCLUSIONS 

This paper presents application of a new weighted particle scheme for conservation laws to the equations of ideal hydro- 
dynamics and magnetohydrodynamics. This scheme has no free parameters which control the physics of the inter-particle 
interaction. There only free-parameter in our scheme is the average number of neighbours that a particle interacts with, and 
this depends only on the number of spatial dimensions. The interaction between particles is entirely described by the source 
terms and fiuxes. The latter are given by the solution of the associated Riemann problem, which correctly treats dissipative 
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Figure 21. Density structure of the disk in XY-plan passing through the centre. The right and left panels show non-magnetised and 
magnetised cases respectively, at t = 100 (top) and at t = 150 (bottom). The density structure in magnetised disk at t = 50 has same 
form as non-magnetised, and therefore is not shown here. 



processes and discontinuous solution without explicit use of artificial viscosity or resistivity. Due to use of Riemann solvers 
and reconstruction methods, our weighted particle method is expected to be more dissipative compared to pure Lagrangian 
SPH without any explicit diffusion terms. However, further quantitative comparison to SPH is required to verify this claim 
in realistic astrophysical problems, where dissipative processes, albeit locally, must be used. 

In our weighted particle scheme, the smoothing length is a property of the particle distribution only, and not the underlying 
solution. As a result, high resolution is obtained in regions with high particle density, which does not need to coincide with high 
mass density regions. This therefore permits similar resolution in both low and high density regions, if the scheme is combined 
with particle refinement methods. In our scheme, the physical meaning of the particle as a fiuid element is lost, and the 
particles should be considered as interpolation points only. Even without the refinement, the mass of the particle can change 
in the course of simulations, though these changes are small in smooth fiows. The advection of a scalar field in our scheme is 
not as trivial as it is in SPH. Namely, for every scalar, a transport equations must be solved which further increase, albeit 
little, both memory and performance footprint of the simulation. If one needs to follow multiple fiuid composition, the matter 



is a bit more complicated, since one will be required to use a consistent multi-fluid advection for chemical species ( Plewa fc 
|Miiller|[l999| ), which in turn further complicates the scheme compared to SPH. Nevertheless, we feel that the advantages of 
the scheme outweigh these disadvantages. 

Our scheme, in principle, permits adaptive particle splitting to increase resolution in the desired regions. However, if not 
done carefully, this may break the regularity of the particle distribution, and therefore introduce substantial errors in flux 
divergence. The magnitude of these errors and their impact on the outcome of the simulation are hard to estimate. The MHD 
test problems with initially non-uniform but properly relaxed particle distribution (j ^5.4|5.5] and |5.6| showed excellent results 
with initially nested and relaxed particle distribution. In principle, if a group of particles is added in the course of simulation 
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Figure 22. The left panel show the amplitude of toroidal field component, and the right panel the magnitude of Maxwell stress, BrBg 
at t = 50 (top) and t = 100 (bottom). The field and the stresses at t = 150 have the structures similar to those at 4 = 100, and therefore 
are not shown here. 



and their neighbours are adjusted to maintain regular distribution, the noise should be small. However, further research is 
required to discover optimal way for particle refinement. 

We showed that the application of the weighted meshless scheme to the equations of ideal hydrodynamics is straightfor- 
ward. We also expect that such scheme is computationally slightly more expensive than SPH. In its optimal form, it requires 
two loops over neighbours, compared to one in SPH: calculation and limit of the gradients of primitive fluid variables, and 
the interaction part. This is in addition to Eq.( |15[ ), which, similarly to SPH, is solved iteratively. The interaction part of the 
scheme requires solution of a Riemann problem between a particle and its neighbours; on average, 32 Riemann problems are 
solved for each particle in 3D. Nevertheless, both HLL or HLLC solvers are only moderately expensive compared to calculation 
of artificial viscosity and conductivity in SPH. However, our scheme has higher memory footprint due to storage of gradients 
in the memory. Nevertheless, combined with the lower number of neighbours usually uses in 3D SPH simulations and a large 
step size in the vicinity of strong shock, which in SPH is limited by artificial viscosity, the performance impact is minimal. 

The strength of the weighted particle scheme becomes clear with its successful application to the equations of the ideal 
MHD. Here, the main problem is the maintenance of V • B = constraint. While it is not clear how to maintain this to 
machine accuracy, if possible at all, in a meshless scheme, this work demonstrated that it is possible to keep the divergence 
under the control by applying hyperbolic-parabolic divergence cleaning method designed for Godunov MHD scheme. This 
MHD formulation is not limited to mess-less schemes, but can also be applied to unstructured grid or moving mesh schemes 
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(Gaburov & Levin 2010, in preparation), in which constraint transport discretisation of the induction equation might be 
difficult or impossible to formulate. 

Finally, we report that in our implementation, the science rate of this meshless MHD scheme is roug hly 10'' particle/s 
in 3D and 2 • 10'* particles/s in 2D on a single 2.7 GHz Core i7 processor core. Being a particle-based scheme, it can also 
be implemented in OpenCL to allow efficient execution on many-core chips, such as GPUs. The science rate of a GPU code 
which we developed and used for the simulations in this paper is 10^ particles/s in 3D, and twice that in 2D on GT200 chip. 
These are the lower values than we expected, and with further tuning and optimisation of the code these rate can be, at least, 
doubled. 



APPENDIX A: PIECEWISE PARABOLIC RECONSTRUCTION 

Approximation of q{x) in the neighbourhood of a particle i is given by a second-order Taylor expansion from the point Xi 

q{x) ^ + {x- Xi)°'q°' + (x - Xi)°'{x - Xi)'^ q^^ . (Al) 

To complete the reconstruction, the coefficients qf and q"^ must be determined. Since the number of the neighbouring particles 
is larger than the number of unknown parameters, this can be done in the least-square sense by minimising the following 
functional ( |Maron fc Howes|2003" l 

Here, Sqij — qj — qi, — x" — xf and Wj is the weight of a particle j, which can be to Wj = w{xj) although other choices 
are possible. The conditions dLi/dqf" = and dd/dq"^ = result in the following set of equations 

= afsr + (A3) 

Qf = q^Sf- + q^'^Sf^-' . (A4) 

Here, Q"''''' = and g"'''^"'^' — u)j^g^fj[^^^.]]^,''j], where the expressions within [..] can be omitted to obtain 

Q", 5°'^ and S""". This system can be solved using methods of linear-algebra, by noticing that both (qf , qf) and (Qf , Qf) 
coefficients can be combined into 9-dimensional vector, and S coefficients into 9x9 matrix. Finally, a parabolic reconstruction 
of an i-particle state at Xij is 

qii;i = qi + Ti \{xij — Xi)"qi + [xij — Xi)°'{xij — Xi)^qf'^ . (A5) 



Here, is a limiter function defined in the same way as for the linear reconstruction (j ]2.2[ ), with one exception: = 0, if an 
extrema occurs within half- vector connecting particles i and j. With such restriction the reconstruction is guaranteed to be 
monotonic. 

This parabolic reconstruction is computationally more expensive compared to the linear one due to large amount of the 
storage and the number of operation required to compute Q and S coefficients, as well as to invert 9x9 matrix. While there is 
advantage of using parabolic reconstruction for Eulerian calculation, at appears to result in little improvement when particles 
move with fluid velocity. 



APPENDIX B: HLL-TYPE MHD RIEMANN SOLVERS 

The one-dimensional MHD equations have the following conservative form 

dt dx ' 

where G = J- — aJJ is a flux in moving frame, W is a fluid state in conservative variables and is a flux in lab-frame: 
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Here, Pt — pth + B^ /2 is the sum of the thermal and magnetic pressures, e — pv^ /2 + eth + B^ /2 is the total energy density, 
and Bx is a constant fleld, implied by V • B = constraint in ID. The initial conditions to this Riemann problem are given 
by specifying left, Ul, and right, Ur, states at the interface at zero time. The flux Q at any time through and normal to 
the interface is given by the Riemann solver. Below, I provide formulae for two commonly used Riemann solvers: HLL and 





HLLD. Both Riemann solvers reduces to the hydrodynamic case if the magnetic field strength is zero. Namely, HLL reduces 
to hydrodynamic variant of HLL solver, and HLLD reduces to HLLC presented in §4.1| 



Bl HLL Riemann solver 



The HLL Rieman solver use a single state to approximate intermediate wave structure ( Harten et al. [1983 1. As a result, it is 
a rather diffusive solver which defuses contact discontinuity even at rest. Nevertheless, it is computationally inexpensive and 
robust Riemann solver, which can be used in pathological cases where other solvers fail. 



The wave structure of the HLL solver is shown in Fig.(Bl I, and the fluxes are given by the following expression 



aJJ* , 



ax < Sl, 

Sl < < Sr, 

Sr < Sr, 



where Tl ~ T{Ul), J^r = J^{Ur) and 



SrJ-l — SlJ'r + SrSl{Ur — Ul) 
Sr — Sl 



is HLL-fiux is the rest frame. The intermediate state is given by the formula 

SrIAr — SlUl — J'R. + J'L 



U = 



Sr — Sl 



(B3) 



(B4) 



(B5) 



The wave speeds are Sl = vcun{vxL,VxR) — Cs and Sr = mSiK{vxL,VxR) + Cs, where Cs = max(cai,, Csij) is the maximal signal 
of the left or right states. More accurate estimates based on Roe-averages are able to reduce diffusion, but the degree of this 
reduction is rather small. If diffusion needs to be minimised, one should rather use intrinsically less-diffusive solver, such as 
linearised Roe-solver, or HLLD Riemann solver described next. 
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B2 HLLD Riemann solver 



From the wave-structure of HLLD solver (Fig B2l, it is clear that contact, Alfven and fast magnetosonic waves are resolved. 
The fluxes at the interface for this solver are given with the following formulae 



^HLLD , 



J'L — axUh, (ix < Sl, 

J^L + (Sl - a^.)Ul - SlUl, Sl < < SI, 

+ {SI — ax)U'l* — [SI - Sl)Kl — SlUl, si < ax < Sm, 

J'R + {Sr — ax)Uji — (Sjj — Sr,)Ui! — SrUr, Sm < ax < Sr, 

J^R + (Sr - ax)UR - SrUr, Sr < a^ < Sr, 

J'R — axUR, Sr < ax- 



(B6) 



The wave speeds are Sl = min(iii,i,, Vxr) — Cs and Sr = max.{vxL,VxR) + Cs, where Cs = max(csL, Csr) is the maximal signal 
speed in left or right state. The speed of the middle wave is defined by 

Pt,r — Pt.l + PVxl{Sl ~ Vxl) — PVxr{Sr — Vxr) 



Sm 



Pl{Sl - Vxl) - Pr{Sr - V^r) 



(B7) 



v*}i — v*if = Sm for K — L, R, and Sl ~ Sm — \Bx\/ \/pl ^-nd Sr — Sm — \Bx\/ ^/'pr- The density of both and ★★-states 
are 

Sk — VxK 



Pk ^ Pk = Ph 



for K = L,R. The next four ★-states are 
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where in [yz] means either y- or z-component of the corresponding 3D field. Following Miyoshi & Kusano ( 2005[ ), if the last 
terms on the RHS result in 0/0 uncertainty, v*y^jf^ = V[yz]K and -B*,,^]^ ~ 0, since in this case there is no shock across Sk- 
Finally, the remaining two ★-states are 



{Sk - VxK)eK - Ptk'"xK -f P}j^Sm + Bx(^k ■ Ba' - ■ B^ 



Sk — Sm 

pi, _ {Sr — Vxr)prPTl — {Sl — Vxl)plPtii + PlPr{Sr — Vxr){Sl — Vxl){vxR — Vxl) 
^ {Sr — ur)pr — {Sl — ul)pl 

and Py*, = P^^ — P^, for K = L , R. The ★★-states are 

** _ VpL'"lyz]L + VPR-"lyz]R + iB[y:,]R ' B[;^]£)sign(_B^) 



^pI + Vpr 



Bu 



X/^Bly^lR + y^B[y^]L + \/P*Lp*R{nyz]R - V[yz]L)sign{Bx) 

Sk ^ BkT \/pI{^k ■ Ba - v** ■ B**)sign(B:r). 
In the last equations, the — and -|- on the right hand side correspond to K = L and R respectively. Finally, v 
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(B13) 

(B14) 

(B15) 
and 



This completes the description of HLLD Riemann solver. It is clear, this solver is computationally more expensive than 
HLL solver due to its ability to also resolve contact and Alfven waves, which is important in many MHD problems. In practice, 
it is proven to be a robust solver for many problems, and tests demonstrate it is at least as accurate as linearised- Roe solver. 
Furthermore, HLLD solver is possible to use with equation of state other than of ideal gas, since the equation of state in 
solver is used implicitly by specifying both the gas pressure and total energy. 
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